** -----------------------------------------------------------------------------
** Calculate p-val for differences in sorting across cohorts
** adjusting for multiple testing using the stepdown method
** -----------------------------------------------------------------------------



** Routine to calculate p-val
** -----------------------------------------------------------------------------

	cap program drop get_pval
	program define get_pval
	
	{

		qui preserve
		local mk = $mk
		di `mk'
		forvalues n=1/$nstat {
			
			qui gsort rep -t_dI
			qui cap drop n
			qui by rep: gen n=_n
			
			local ix = index[1]                                  // reference of next max tstat
			global tdI`mk'_`ix' = t_dI[1]                            // next max tstat
			global dI`mk'_`ix'  = dI[1]                              // corresponding difference
			global se`mk'_`ix'  = abs(dI[1])/t_dI[1]                 // and its se
			qui count if ${tdI`mk'_`ix'}<t_dI & rep>0 & n==1         // prob larger tstat for adjusted pval
			global ap`mk'_`ix' = r(N)/$nrep                           // adjusted p-val
									
			qui count if ${tdI`mk'_`ix'}<t_dI & rep>0 & index==`ix'  // prob larger tstat for unadjusted pval
			global up`mk'_`ix' = r(N)/$nrep                           // unadjusted p-val
			
			qui drop if index==`ix'

		}
		qui restore
		
	}

	end	


	
	
** Test for changes in assortativeness across cohorts
** -----------------------------------------------------------------------------

{
	
	local l1 "30"
	local l2 "40"
	local l3 "50"
	local l4 "60"
	local l5 "70"

	local c2=5			// 70s cohort
	forvalues c1=1/4 {	// earlier cohorts

		qui use "$logfolder\bsTests_`l`c2''vs`l`c1''.dta", replace
		global c1 = 1900+10*(`c1'+2)
		global c2 = 1900+10*(`c2'+2)

		qui sum rep
		global nrep = r(max)

		** p-values, all 30 indices together
			global nstat = 30
			global mk    = 0
			qui get_pval

		** p-values, each market separately
			global nstat = 7	// number of joint tests
			forvalues mk=1/4 {	// markets
				qui preserve
				qui keep if inrange(index, (`mk'-1)*7+1, `mk'*7)
				global mk=`mk'
				qui get_pval
				qui restore
			}
			
			global nstat = 2	// number of joint tests
			global mk    = 5	// market
			preserve
			qui keep if inrange(index, 29, 30)
			qui get_pval
			restore
		
		** print estimates
		local c=10*(`c1'+2)
		
		di _newline
		di "------------------------------------------------------------------------"
		di "Cohort ${c2} vs `c'"
		di "Columns: C&PG vs SC,  PG vs C, C vs SC, C&PG vs others, 5 educ levels" _newline
		di "Odds ratio          ",     %6.3f ${dI0_1}    , "   ",     %6.3f ${dI0_8}     , "   ",     %6.3f ${dI0_15}    , "   ",     %6.3f ${dI0_22}    
		di "                    ", "(" %3.2f ${ap0_1} ")", "   ", "(" %3.2f ${ap0_8}  ")", "   ", "(" %3.2f ${ap0_15} ")", "   ", "(" %3.2f ${ap0_22} ")"
		di "                    ", "[" %3.2f ${ap1_1} "]", "   ", "[" %3.2f ${ap2_8}  "]", "   ", "[" %3.2f ${ap3_15} "]", "   ", "[" %3.2f ${ap4_22} "]" _newline
		di "L1 index            ",     %6.3f ${dI0_2}    , "   ",     %6.3f ${dI0_9}     , "   ",     %6.3f ${dI0_16}    , "   ",     %6.3f ${dI0_23}    
		di "                    ", "(" %3.2f ${ap0_2} ")", "   ", "(" %3.2f ${ap0_9}  ")", "   ", "(" %3.2f ${ap0_16} ")", "   ", "(" %3.2f ${ap0_23} ")"
		di "                    ", "[" %3.2f ${ap1_2} "]", "   ", "[" %3.2f ${ap2_9}  "]", "   ", "[" %3.2f ${ap3_16} "]", "   ", "[" %3.2f ${ap4_23} "]" _newline
		di "L2 index            ",     %6.3f ${dI0_3}    , "   ",     %6.3f ${dI0_10}    , "   ",     %6.3f ${dI0_17}    , "   ",     %6.3f ${dI0_24}    
		di "                    ", "(" %3.2f ${ap0_3} ")", "   ", "(" %3.2f ${ap0_10} ")", "   ", "(" %3.2f ${ap0_17} ")", "   ", "(" %3.2f ${ap0_24} ")"
		di "                    ", "[" %3.2f ${ap1_3} "]", "   ", "[" %3.2f ${ap2_10} "]", "   ", "[" %3.2f ${ap3_17} "]", "   ", "[" %3.2f ${ap4_24} "]" _newline
		di "Weighted L ratio    ",     %6.3f ${dI0_4}    , "   ",     %6.3f ${dI0_11}    , "   ",     %6.3f ${dI0_18}    , "   ",     %6.3f ${dI0_25}    , "   ",     %6.3f ${dI0_29}    
		di "                    ", "(" %3.2f ${ap0_4} ")", "   ", "(" %3.2f ${ap0_11} ")", "   ", "(" %3.2f ${ap0_18} ")", "   ", "(" %3.2f ${ap0_25} ")", "   ", "(" %3.2f ${ap0_29} ")"
		di "                    ", "[" %3.2f ${ap1_4} "]", "   ", "[" %3.2f ${ap2_11} "]", "   ", "[" %3.2f ${ap3_18} "]", "   ", "[" %3.2f ${ap4_25} "]", "   ", "[" %3.2f ${ap5_29} "]" _newline
		di "Normalized trace    ",     %6.3f ${dI0_5}    , "   ",     %6.3f ${dI0_12}    , "   ",     %6.3f ${dI0_19}    , "   ",     %6.3f ${dI0_26}    , "   ",     %6.3f ${dI0_30}    
		di "                    ", "(" %3.2f ${ap0_5} ")", "   ", "(" %3.2f ${ap0_12} ")", "   ", "(" %3.2f ${ap0_19} ")", "   ", "(" %3.2f ${ap0_26} ")", "   ", "(" %3.2f ${ap0_30} ")"
		di "                    ", "[" %3.2f ${ap1_5} "]", "   ", "[" %3.2f ${ap2_12} "]", "   ", "[" %3.2f ${ap3_19} "]", "   ", "[" %3.2f ${ap4_26} "]", "   ", "[" %3.2f ${ap5_30} "]" _newline
		di "$\chi^2$            ",     %6.3f ${dI0_6}    , "   ",     %6.3f ${dI0_13}    , "   ",     %6.3f ${dI0_20}    , "   ",     %6.3f ${dI0_27}    
		di "                    ", "(" %3.2f ${ap0_6} ")", "   ", "(" %3.2f ${ap0_13} ")", "   ", "(" %3.2f ${ap0_20} ")", "   ", "(" %3.2f ${ap0_27} ")"
		di "                    ", "[" %3.2f ${ap1_6} "]", "   ", "[" %3.2f ${ap2_13} "]", "   ", "[" %3.2f ${ap3_20} "]", "   ", "[" %3.2f ${ap4_27} "]" _newline
		di "Minimum distance    ",     %6.3f ${dI0_7}    , "   ",     %6.3f ${dI0_14}    , "   ",     %6.3f ${dI0_21}    , "   ",     %6.3f ${dI0_28}    
		di "                    ", "(" %3.2f ${ap0_7} ")", "   ", "(" %3.2f ${ap0_14} ")", "   ", "(" %3.2f ${ap0_21} ")", "   ", "(" %3.2f ${ap0_28} ")"
		di "                    ", "[" %3.2f ${ap1_7} "]", "   ", "[" %3.2f ${ap2_14} "]", "   ", "[" %3.2f ${ap3_21} "]", "   ", "[" %3.2f ${ap4_28} "]"
		di "------------------------------------------------------------------------"
		
	}

}


